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Black hole mergers: the first light 
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ABSTRACT 

The coalescence of supermassive black hole binaries occurs via the emission of gravitational 
waves, that can impart a substantial recoil to the merged black hole. We consider the energy 
dissipation, that results if the recoiling black hole is surrounded by a thin circumbinary disc. 
Our results differ significantly from those of previous investigations. We show analytically 
that the dominant source of energy is often potential energy, released as gas in the outer disc 
attempts to circularize at smaller radii. Thus, dimensional estimates, that include only the 
kinetic energy gained by the disc gas, underestimate the real energy loss. This underestimate 
can exceed an order of magnitude, if the recoil is directed close to the disc plane. We use three 
dimensional Smooth Particle Hydrodynamics (SPH) simulations and two dimensional finite 
difference simulations to verify our analytic estimates. We also compute the bolometric light 
curve, which is found to vary strongly depending upon the kick angle. A prompt emission 
signature due to this mechanism may be observable for low mass (10^ Mq) black holes whose 
recoil velocities exceed ~ 10^ km s"'. Emission at earlier times can mainly result from the 
response of the disc to the loss of mass, as the black holes merge. We derive analytically the 
condition for this to happen. 

Key words: black hole physics — accretion, accretion discs — hydrodynamics 



1 INTRODUCTION 

Supermassive black hole binaries are predicted by hierarchical 
galaxy formation models and are a likely consequence of observed 
galaxy mergers. Howeve r, only a h andful of these binaries have 
been directly observed ( Rodriguez et al. 2006) and their dynami- 
cal evolution is still uncertain. If binaries coalesce on a time-scale 
shorter than the age of the universe, mergers can be an important 
ingredient in the evolution and growth of supermassive black holes. 
Mergers also emit low frequency gravitational waves whose detec- 
tion is one of the prime goals of proposed experiments such as the 
Laser Interferometer Space Antenna (LISA). What remains uncer- 
tain is whether detectable electromagnetic emission occurs either 
prior to or in the immediate aftermath of the coalescence. 

Observations and numerical simulations strongly suggest that 
mergers of gas-rich galaxies result in the inflow of gas into the nu- 
clear region, where it is likely to co-exist with a newly formed black 
hole binary (' Escala et alTllfoOSi : lOotti et alj|2007l : (Callegari et al.l 
|2009). The gas may have an active role and mediate and speed u| 



2UUyj. ine gas may nave an active role ana mediate ana speed up 
the coalescence i Ivanov et al. 19991: Armitage & Nataraiaiill2002l : 



ICuadra et alj|2009l : iLodato et al.ll2009h ~ However - and irrespec- 
tive of whether the gas is dynamically important - if gas is present 
then perturbations to the gas during the merger may give rise to 



observable electromagnetic signals that can either precede or fol- 
low the gravitational wave signal. Observations of these signals can 
give us evidence for the existence of mergers, improve upon LISA'S 
limited ability to spatially localise sources, and teach us about the 
astrophysical phenomena involved in the process. 

General relativity predicts that gravitational waves emitted by 
the shrinking binary during the final coalescence carry away a non- 
zero net linear mom entum, so that the centre of mas s of the merged 
black holes recoils ( |Pereslll962l : iMiller et al]|2007l) . A number of 
authors have used semi-analytic methods or N-body simulations 
to estimate the influence of the recoil on a surrounding thin disc 
jLip yai et alj200 3: IShields & Bonning|l2008l : ISchnittman & KrolikI 
2008). Hydrodynamic simulations for a thick disc - a configura- 
tion that we do not co nsider here - have also been performed by 
iMegevand etal.l ( l2009l) . 

Here we revisit the analytic problem (emphasizing the im- 
portance of different physical effects from those discussed pre- 
viously), and perform numerical simulations to investigate how 
a thin, non self-gravitating disc reacts to the birth of a central 
black hole via merger. The initi al disc geometry that we con sider 
is based upon that discussed bylArmitage & Nataraianl ( |2002|) and 
iMilosavlievic & PhinnevI ( l2005h . We assume that, at a time signif- 
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icantly prior to the final coalescence (when the binary semi-major 
axis was > 10^ Schwarzschild radii), the binary was surrounded by 
a geometrically thin disc whose plane coincided with that of the 
binary. While the binary remains in tidal contact with this disc, en- 
ergy and angular momentum can be exchanged between the two. 
This excha nge provides an additional source of heat for the inner 
disc (lLod"a to et al. 2009), and alters its long term viscous evolution 
( |Pringlelll991D . Immediately prior to the coalescence, however, the 
rapid decay of the binary orbit due to gravitational wave emission 
leads to a decoupling of the binary from the inner disc. As a re- 
sult, it is a reasonable first approximation to model the eff'ect of a 
kick on an initially unperturbed, axisymmetric circumbinary disc. 
We also ignore any effects associated with gas that has survived in 
circumprimary or circumsecondary discs up until near the moment 
of co alescence. Such gas, if present, may produce observable signa- 
tures ( lArmitage & Nataraiaj20o3 : lLodato et alj2009l : lchang et all 
l2009l) . but it is decoupled dynamically from the gas in the circumbi- 
nary disc. The effect of the recoil on the circumbinary disc can 
therefore be considered independently from the fate of gas bound 
to either of the individual black holes. 

In this paper, our main focus is on the effect of the black hole 
recoil rather than the a lmost instantaneous mass loss t hat also ac- 
companies the merger jMilosavlievic & Phirmevll2005h . We study 
the consequence of velocity perturbations in both the linear (kick 
velocity smaller than circular velocity) and non linear regimes. We 
calculate the magnitude of energy dissipation and assess the depen- 
dence of this potentially observable quantity on the recoil geometry. 
In particular, we focus on off-plane kicks, which are expected to be 
the most common. 

The paper is organised as follows. In Section|2]we analytically 
investigate the properties of the disc after the kick and the dissipa- 
tion of the extra energy in the innermost region. These estimates 
serve as guidance for our numerical simulations, described in Sec- 
tion [3] In Section |4] we present our numerical results, which we 
discuss in Section|5] Finally, in Section|6]we draw our conclusions. 



2 ANALYTIC ESTIMATES 

2.1 Disc topology: bound and unbound regions 

We first consider the reaction to the kick of particles at radius R. 
Initially, the particles rotate anti-clockwise around the black hole 
in circular orbits with Keplerian velocity = yJGM/R (j>, where 
(f) is the angle between the positive i/-axis and the particle radius R, 
M is the mass of the merged black hole and G is the gravitational 
constant (see Fig.[T]l. The angular velocity is £1 = -\JgmJW. 

In the frame of the moving central object the recoil results in 
each particle receiving an an additional velocity V. The kick veloc- 
ity direction makes an angle 6 with the disc plane. For simplicity, 
we assume that the projection onto the disc plane, Vcos 9, is in the 
direction of the positive x-axis (see Fig. [T). In cylindrical coordi- 
nates the three components of the initial particle velocity are then 
Vr = - V cos 6 sin 0, = - V cos 9 cos (p and = V sin d. 

If the disc extends to sufficiently large radii, there is a radius 
R, at which V = V^. In the following, it is convenient to adopt the 
dimensionless radial coordinate r = R/Ry = (V/V]^)'. The specific 
energy of a particle then reads 



-i^^(l+2V^c 



(1) 



A particle can escape from the system (e > 0) or it can remain 
bound to it (6 < 0), depending on its distance from the hole and on 
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Figure 1. Face-on view of the disc showing the the Keplerian velocity 
and the projection of the kick velocity onto the disc plane V cos 6. The an- 
gle between the positive y-axis and the particle's position is 0. The closed 
solid lines mark the border between the bound and unbound regions for 
kick angles between 0° (a kick directed into the disc plane) and 90° . As 6 
approaches 90 degrees, the bound region becomes increasingly azimuthally 
symmetric. 



the direction of the kick velocity relative to its original Keplerian 
velocity. In particular at small radii. 



r < rt, = (- cos6 + VcosS- + l] , 

all particles within a given annulus are bound, while for 



(cos 9+ Vi 



cos 9^ + 



(2) 



(3) 



all particles are unbound. The dimensionless radii = Rb/Ry and 
''ub = Riib/Ry correspond to the actual radii R\, and R^b, respectively. 
In the radial range between rj, and r^^b, particles at radius r are bound 



only within a limited azimuthal range 
k), where 

r- I 



'b < ' 



, (with ^ . 



(4) 



2 Vrcos 6 

In Fig. [T] we show the boundaries of the bound region for different 
kick angles. As 9 increases, the bound region becomes increasingly 
azimuthally symmetric, while the radial extent of the bound region 
between rb and r^b decreases. For a perpendicular kick rb = r\,b = 1 • 
In our analytic estimates we consider only the portion of the 
disc that is still bound after the kick, since matter that is unbound 
by the recoil simply leaves the system and it is not energetically 
important. 

2.2 Angular momentum 

In the bound portion of the disc (where e < 0), particles are excited 
on elliptical orbits with eccentricity e- = 1 -I- 2J-e/(GM)-, and with 
specific angular momentum 



l.(l-V7c 



(5) 



where /i; = RVt^ is the Keplerian specific angular momentum for 
that radius. In general j t Jk, and as a consequence rearrangement 
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Figure 2. The mean circularisation radius (equation|7) in units of i?b (equa- 
tion |2) as a function of the initial location R in the same units. The black 
solid thin line {R = R^) divides the region where outward spreading occurs 
(upper part R^ > R) from the region where accretion occurs (R^ < R). 



of gas must occur after the recoil as the gas seeks a new minimum 
energy configuration. If a particle conserves its angular momentum 
then it must circularize at a radius r^.p = R^.p/Ry given by 



(l- -fr cos 6 cos 0) + r sin^ Q. 



(6) 



To understand the behaviour of the disc after the kick, we calculate 
the mean circularization radius for a ring, = R^/Ry, 



r 20b J-^b \ Jk 



rcos' 6 



sin Icpb 
20b 



(1 + rsin" 6) + 
sin 0b 



(7) 



where the integration limits are such that 0b = 't. for r < r^,. Fig.|2] 
shows that while matter R < R\, circularizes close to its orig- 
inal location for any kick angle, R^ depends strongly on the kick 
angle for R > R^,. For high latitude kicks (6 t 50°), the behaviour 
outside and inside R^ is similar. For smaller inclination angles an 
increasing amount of matter circularizes at smaller radii, and only 
a small fraction of matter expands toward larger radii. This implies 
that there must be a net release of potential energy. 

For out-of-the plane kicks, the recoiling matter acquires in 
general an azimuthal component in its new angular momentum 
(equation |5]l. Its magnitude varies with radius, so we expect the 
disc to be warped initially. After the extra energy associated with 
the warp is dissipated, the disc must lie in a plane, however tilted 
with respect to its original orientation. 



2.3 Energy dissipation 

In the previous section, we showed that particles that remain bound 
after the kick do not, in general, have the same specific angular 
momentum as they had prior to the kick. Initially their orbits are 
eccentric. In a fluid disc, orbit crossing leads to energy dissipation 
and circularization of the disc into a new circular equilibrium con- 
fiauration. 



We estimate the total energy available for dissipation by as- 
suming that each particle that remains bound circularizes at the ra- 
dius appropriate to its post-kick angular momentum. We write this 
as, 



(8) 



where E is the initial post-kick energy of the bound disc (including 
the gain or loss of kinetic energy resulting from the kick) and E^ 
is its energy after circularization. The initial energy is obtained by 
integrating the specific energy (equation [T) over the bound region 
of the disc. 



(El,RdRd(/). 



(9) 



Here E is the disc surface density, Ri„ is the inner radius of the disc, 
and it is to be understood that the limits on the angular integral are 
such that 0b = TT when R si R\,. 

After circularization the total energy of the gas within the disc 

is 



E, 



(10) 



where (r/r^^p) is given by equation The integral over is ana- 
lytic and can readily be evaluated using a symbolic algebra pack- 
age. However the simplest general form we have found is extremely 
lengthy, and we do not write it here. 

To assess which parts of the disc contribute the most to E^iss 
we now specialize to a disc with a power-law surface density profile 
E = Y,\r^'', where Ev = E(r =1). Dimensional argument^ suggest 
that the change in kinetic energy is of the order of (V-/2)Sv/?y. 
Hence we write the differential contribution to the total energy re- 
lease as. 



d£diss 
dr 

where, 

g{r) -- 

I{r) : 



^\ R^- 



20b r-<' 

20b J-(6b Uk 



r-2V^cos0^^-I+J(r) 



0b 



(11) 

(12) 
(13) 



The function g(r) is a measure of how well the simple estimate 
(l/2)V-Si/^v captures the expected total energy release. This func- 
tion is plotted for different values of the kick angle (and p = 1.5) in 
Fig. [3] As 6 decreases, the behaviour around r = I is increasingly 
dominated by I and it develops a sharp peak at r = 1. For 6 = 0°, 
the right-hand side of equation jilt diverges as \r- 1|"'. For d t 0° 
we integrate over radius to obtain the final estimate of iidiss- The 
dependence of E^jiss on kick angle is shown in Fig.|4] We find that 
regardless of the surface density distribution, E^i^s is almost 3 or- 
ders of magnitude larger for a kick close to the disc plane than for 
a perpendicular kick. Moreover the assumption made in previous 
papers that (l/2)V-'I.vRy would fairly measure E^i^s proves to be a 
substantial underestimate of the energy released for all kick angles 
e < 50°. 



For example, ISchnittman & KrolikI j2008l) estimate the total energy re- 



lease to be 4{V^/2)I.y Ri. 
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Figure 3. Differential energy to be dissipated is plotted as a function of 
radius (equation II U . The energy is expressed in units of (\ l2)V^'LvR^, 
while the radius is in units of R\. Different curves correspond to different 
kick angles as labelled. The surface density profile is Z oc r"'/^. 
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Figure 4. Total energy available to be released (equation [8) is plotted as a 
function of the kick angle for different power-law indexes for the surface 
density distribution. The energy is expressed in units of (l/2)y^ZvRy. 



2.4 Contributions to E^^s 

The energy available as the disc adjusts to a new equilibrium con- 
figuration can be considered to have two sources. First, there is 
an immediate change to the kinetic energy of particles because of 
the recoil. When integrated over the bound region the total energy 
change from this effect is positive for > 4°. Second, there is the 
release of potential energy that is liberated when gas in the disc 
loses angular momentum and circularizes at smaller radii. In the 
remainder of the paper we will refer to the release of potential en- 
ergy on a short time scale as "accretion energy". This should not be 
confused with the (much slower) energy release due to viscous disc 
accretion (described as "accretion energy" by Schnittman & Krolik 
2008), which we do not model. 

To quantify how the relative importance of the two energy 
sources varies with radius and with kick angle we rewrite the post- 
kick specific energy as 



- (2 V'' cos 6 cos (f>- rj 
and define the prompt change to the energy of the disc as. 



(14) 



(15) 



We note immediately that if (? = 90° then e^ick is simply equal 
to V- 12. Integration over the bound region of the disc then yields 
£kick = ( 1 /2)Mdisc V^, where Mdisc is the total mass within the bound 
region. 

For arbitrary kick angles one can readily show that the contri- 
bution to iiuck for < is given by the same expression, 

A£kick = ^AMdi^cV", (16) 

but where AMdisc is now the mass within the radius where the entire 
annulus remains bound to the hole. At larger radii, however, this is 
no longer true. For R\, < R < R^b the differential contribution to the 
prompt energy change dEkick/dR can be either positive or negative, 
depending upon the radius and upon the specific kick angle. This 
is illustrated in Fig. |5] for the case of a kick at 6 = 30°. When the 
integration is extended across the entirety of the bound disc the 
result is that Euck * (l/2)Mi,,^V^ for # 90°. 

Finally we can compare the relative importance of the prompt 
kinetic energy to the part of £diss that arises from potential energy 
release. In general, we find that for R > R], the gain in kinetic energy 
is either negative or (when positive) negligible compared to the ac- 
cretion energy, whereas for R ^ Rt it can be substantial (Fig.|5]l. As 
we will discuss subsequently the physical size of R^ is often quite 
large, so for a small disc that does not extend out to Rf, it may be 
a reasonable approximation to assume that the disc receives an in- 
jection of kinetic energy that is subsequently dissipated. For a large 
disc, on the other hand, any kick that is not nearly perpendicular 
to the disc plane produces a wholesale rearrangement of the gas. 
The total energy dissipation in this regime is dominated by accre- 
tion energy, and detailed hydrodynamic simulations are needed to 
assess accurately the magnitude and time scale of energy release. 

When we consider the energies integrated over the whole 
bound disc, we find that as 6 decreases the fraction of isjiss that 
is due to release of potential energy increases. This is because the 
recoil creates larger gradients in the specific energy and angular 
momentum of the particles. Two effects result. First, the region be- 
tween fb and fub grows and, second, within that region more par- 
ticles have a new lower specific angular momentum (see Fig. 
Therefore, hydrodynamic simulations are needed, especially, for 
kicks grazing the plane and above all, for in-plane kicks, where 
the analytic estimate of Ediss fails. 



2.5 Effect of mass variation of the post-remnant black hole 

Not only momentum, but also energy is carried away by the grav- 
itational waves that cause the final merger of a black hole bi- 
nary. This results in a net mass loss for the system and the fi- 
nal black hole born by merger has a smaller mass than the sum 
of the masses of the two progenitors. This instantaneous mass 
loss of the central attractor has a n effect on the surrounding 
gas (Milosavljevic & Phinnev 2005; Sch nittman & KrolikI |2008| - 
O'Neill et al. 2009; Megevand et al. 2009). We show here that this 
effect deposits substantial energy only within the innermost ~ 10'' 
Schwarzschild radii for V > 100 km s"'. 

The argument is as follows. We first estimate the energy dis- 
sipated in the gas as the disc recovers a circular equilibrium con- 
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Figure 5. Differential energy available for dissipation (solid line) and 
change rn total kinetic energy imparted directly by the kick (dashed line) 
is plotted as a function of radius, for 6 = 30°. As previously, energy is in 
units of (I/2)V-ZvRy> '^e radius is in units of R\, and S oc r^^l^. 
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Figure 6. Radius within which Ae„, dominates over 12 as a function of 
the kick velocity V. The dashed line is an average over the spins, assuming 
isotropic distribution of spin orientation (it amounts to the curve for the 
kick velocities due only to different mass ratios for the BHs). The solid line 
assumes equal magnitude spins, aligned with the orbital angular momentum 
vector but pointing in opposite directions. Finally, the dot-dashed line is for 
equal magnitude spins, in the direction perpendicular to the orbital angular 
momentum and anti-aligned. 



figuration around a lighter central object, neglecting the effect of 
recoil. The effect of the mass loss 6M is to put particles in ec- 
centric orbits, since it instantaneously decreases their binding en- 
ergy. The new specific energy of a ring initially at radius R is 
ei = -GM/(2R)(l - 2SM/M). In the circularization process, the 
particles settle on a circular orbit at a larger radius R^, given by 
angular momentum conservation: i?c = ^/(l - SM/M). At this final 
radius the specific energy is 6f = -GM/{2R){1 - SM/M)^. Thus, 
the induced epicyclic motions dissipate an energy 



A6„, 



(17) 



-vli-] 

We now assess the relative importance at a given radius of the 
effects of the mass loss compared to that of the recoil. Equation l ll7t 
shows that Acm oc R ' is a decreasing function of radius. Since the 
contribution of mass loss turns out to be important only within ^i,, 
we can compare Ae^ with the energy imparted by the recoil a 
We find that there is a radius R,^ < R], within which Aen, > V^/2, 



M I 



(18) 



where R^ = 2GM/c^ is the Schwarzschild radius and c is the speed 
of light. Both 6M/M and V depend on the mass ratio and on the 
magnitudes and directions of the spins of the merging black holes. 
To evaluate R^ we us e the analytic expres s ions, calib rated against 
simula tions, given by iTichv & Marronettil (l2008l) and iMiller et"aD 
respectively. The resulting value of Rm/Rs is shown in 
Fig-Hi for three different configurations of the black hole spins (oi , 
02). For some choices of the parameters, two different magnitudes 
for the mass loss are associated with the same V, explaining the 
double valued behaviour of rJ^ The figure shows that the larger 
the kick, the smaller the radius of influence of the mass loss, ap- 
proaching the innermost stable orbit for the most extreme recoils. 



While SM depends only on the total energy emitted in gravitational 
waves, V depends also on the degree of asymmetry of the emission. There- 
fore, two different configurations at merger can release different total energy 
but carry away the same momentum. 



In our simulations, we verify that for our set-up the dissipation 
from mass loss is indeed negligible at all times. We will show an 
example of how the lightcurve from the recoil compares with that 
of the mass loss for an in-plane kick in Section |4!4l 



3 SIMULATIONS 

Our analytic predictions are based upon a particle rather than a fluid 
model for the disc. Thus, we have run a set of simplified numerical 
hydrodynamic simulations, in order to check these predictions, and 
to obtain the explicit form of the energy dissipation rate within the 
disc as a function of time. We discuss first the special cases of per- 
pendicular (d = 90°) and in-plane (6 = 0°) kicks. The perpendicular 
case is exactly axisymmetric, while the in-plane case is also a two- 
dimensional problem (in r and 0) in the limit of an infinitesimally 
thin disc. Given their symmetry, these cases can be simulated effi- 
ciently using Eulerian finite difference methods. The general case 
of out-of-plane kicks and the 6 = 0° case when the disc has a finite 
thickness are fully three-dimensional, and correspondingly more 
difficult. We treat these cases using t he Lagrangian Smoothed Par- 
ticle Hydro d ynamics (S PH) m ethod l lGingold & MonaghanI [19771: 
lBenzlll990l : iMonagha 3 Il992h . We compare the two numerical 
methods by simulating the 6 = 90° case with both codes. 



3.1 Initial set-up 

We investigate a thin disc, extending from r = 0.1 to r = 10, 
in order to encompass the three regions described in Section [TTI 
The radial distribution of surface density and sound speed, c,, are 
power-laws 



S(r) oc r-" 
cAr) oc r-^'\ 



(19) 
(20) 



with the normalization of the sound speed chosen such that the disc 
aspect ratio H/R = 0.05 at the characteristic radius Ry. The initial 
conditions of the simulations are isothermal with a fixed (in time) 
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but radially dependent sound speed. In the finite difference runs 
the sound speed profile at later times remains fixed (in Eulerian co- 
ordinates), whereas in the SPH run the internal energy per unit mass 
of each particle remains constant (i.e. the sound speed is fixed in 
a Lagrangian sense). In either case the physical assumption is that 
the energy dissipated as the disc readjusts to its new equilibrium 
configuration is radiated away rapidly. Since we do not include the 
effects of the disc self-gravity, the surface density normalization is 
arbitrary. We can freely scale our results to represent a particular 
choice of the disc to black hole mass ratio q = M^i/M. 

The unperturbed — pre-kick — disc velocity field is set to 
give centrifugal equilibrium in the black hole gravitational field (as- 
sumed to be Newtonian). The correction arising from radial pres- 
sure gradients within the disc is taken into account, even if it is 
small for the thin discs that we consider. Working in the frame co- 
moving with the black hole, the initial (perturbed) velocity field is 
a superposition of Keplerian anticlockwise rotation with a kick of 
magnitude V, directed at an angle 6 from the disc. 

3.2 Numerical implementation 

As mentioned above, we performed numerical simulations using 
two very different schemes. We have compared the corresponding 
results in order to check that they are numerically reliable. In this 
section we give an overview of both codes, and in particular note 
those features that are most important for understanding our con- 
clusions. 

We model the e ffect of an in-plane or p erpendicular kick us- 
ing the ZEUS code jStone & Normadll992l) . ZEUS is a Eulerian 
finite-difference code. For in-plane kicks we use cylindrical po- 
lar co-ordinates (r, (p), assuming symmetry in the z direction. For 
perpendicular kicks we use spherical polar co-ordinates (r, 0), as- 
suming symmetry in 0. We impose outflow boundary conditions at 
the inner and outer edges of the grid. Outflow boundary conditions 
are implemented in ZEUS by setting the gradient of all physical 
quantities to zero at the boundary. This is exact only for supersonic 
outflow. Exploratory runs showed that for in-plane kicks enough 
accretion occurs to bring subsonic gas into contact with the inner 
boundary. This leads to unphysical reflections. To ameliorate this, 
we additionally impose a near-vacuum in the innermost active grid 
zone by sweeping away gas above some small threshold on every 
time step. 

SPH is a Lagrangian code to solve the hydrodynamics equa- 
tions. Our version of the code self-consistently includes the so- 
called V/z terms needed to ensure energy conservation (Monaghan 
l2007h . where h is the smoothing length, which sets the effective res- 
olution length of the simulation. The smoothing length is adjusted 
adaptively in the code, ensuring higher resolution in high density 
regions. For a disc with surface density and sound speed profiles as 
in equations l |19t and l l20t . with p = 3/2, the azimuthally averaged 
smoothing length (h) oc H. Therefore, this disc vertical structure 
is equally resolved at each radius. We have run SPH simulations 
using different total number of particles A' in order to check for 
convergence. We have found that our results, and in particular the 
expected luminosity arising from the kick, converge at high reso- 
lution, for N > 4 million. In the following we present only results 
corresponding to our highest resolution simulations, which employ 
either 4 or 8 million particles. 

Both codes use artificial viscosity to capture shocks. In our 
simulations the stress tensor associated with artificial viscosity en- 
ters the momentum equation but is neglected in the energy equa- 
tion (indeed, no energy equation is required in our isothermal runs). 



Nonetheless, we can compute the implied rate at which viscosity is 
dissipating kinetic energy in shocks and compressions. Use this to 
measure the spatial and temporal dependence of the energy dissi- 
pation caused by the kick. 

In ZEUS, in one spatial dimension (say r) the implied rate of 
change of the internal energy per unit volume e is, 



dt ~ ^\ Ar 



(21) 



where Au, = - v^j is the change in velocity across one cell of 
width Ar and. 



Q 



Cp(Av,y ifAv<0 
otherwise. 



(22) 



The dimensionless constant C controls the numerical width of 
shock fronts. We use C = 2. Following the operator-split approach 
used throughout ZEUS we calculate the full heating rate by adding 
independent one-dimensional terms in r and (p (for in-plane kicks) 
or r and 6 (for perpendicular kicks). 

In the SPH code the fundamental quantity is the energy per 
unit mass u. The implied rate of change of u due to artificial vis- 
cosity is 



du Q 
at p 



(23) 



where the divergence of the velocity is calculated over a smoothing 
length h and the artificial viscosity term Q is 



Q 



aCghV ■ V + y 




(V-v)2 ifV-v<0 
otherwise. 



(24) 



We use the viscosity switch described in Morris & Monaghan 
(1997), where the parameter f3 = 2a and a varies from 1 to 0.01 
so that dissipation is enhanced when the particle enters the shock 
region and reduced otherwise. 

We note that the quadratic terms for artificial viscosity are es- 
sentially equivalent in the two codes, although the SPH implemen- 
tation is not operator split. To calculate the total luminosity that is 
released (thereafter labelled "implied dE/dt"), we sum the implied 
rate of energy dissipation (equation d22b or equation i24i ) over vol- 
ume for ZEUS, and over mass for SPH. 



4 NUMERICAL RESULTS 

In the following, results are presented in code units in which the 
unit length is taken to be 



R^j = GM/V^ 

and the unit time is the dynamical time 



fv 



/GM\ 



(25) 



(26) 



I V3 /' 

at Ry, where V is the kick velocity and M the black hole mass. A 
unit time in code units thus translates to f = 155 yr for a black hole 
mass of M = 10'' Mq and a kick velocity V = 300 km s"'. The 
scaling for the energy dissipation rate is, 

dE 



\ ?codc l\G I dt 



(27) 



where dE/dt is the implied dissipation rate in code units, q = 
Mj/M is the actual mass ratio between disc and black hole and 
Icode = 6 X lO"* is a fiducial value that we use when plotting the 
derived energy dissipation rate. For the same kick parameters given 
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Figure 8. Compaiison of the implied energy dissipation rate evaluated from 
ZEUS (solid curves) and SPH (long dashed curves) simulations for perpen- 
dicular kicks. The upper curves show results for steep surface density pro- 
files {p = 3/2), the lower curves results for p = 3/5. The dotted curves show 
the portion of the implied SPH dissipation that arises from the quadratic 
term in the artificial viscosity. 

above and q = ^codc. a typical code luminosity of lO"* then cor- 
responds to a physical luminosity of 3.6 x 10""' erg s"'. Note that 
the peak luminosity depends upon the ratio of the disc mass to the 
black hole mass and ~ very sensitively - on the kick velocity, but 
that there is no dependence upon the black hole mass. 



4.1 Results for perpendicular kicks 

To model perpendicular kicks we use ZEUS in its axisymmetric 
spherical polar mode. For these runs the computational grid extends 
from r = 0.025 to r = 12, and < 6 < n. The disc is set up in 
vertical hydrostatic equilibrium with the mid-plane in the equatorial 
plane. The disc extends from r = 0. 1 to r = 10. All luminosities are 
scaled to represent a disc mass ratio ^code = 6 x 10""*. We use 200 
mesh points logarithmically spaced in r, and 200 mesh points in 9. 

Fig.|7]shows the evolution of the density in the (r, z) plane for a 
disc with a surface density profile S oc r^^^^. There is a clear distinc- 
tion between the behaviour of the gas in the bound and unbound re- 
gions, which in this case are radially separated at r = rb = ;■„(, = 1- 
The bound gas is distorted by the kick into a bowl-like shape, which 
oscillates before settling back into the equatorial plane. The un- 
bound gas departs on an essentially ballistic trajectory. A low den- 
sity streamer of gas that is marginally bound persists out to quite 
late times. 

The energy dissipation rate is plotted as the solid lines in 
Fig. [8] The luminosity peaks at about t/ty = 0.4 for both of the 
surface density profiles that we simulated (p = 3/2 and p - 3/5), 
before declining slowly out to late times. The oscillations seen after 
the peak are present in both the SPH and the ZEUS results and also 
in ZEUS simulations run at different resolutions. Thus they appear 
to be physical. Since they do not appear for other kick angles they 



are probably associated with the dissipation of vertical shear in the 
90° case. 

4.1.1 Numerical tests 

Fig. [8] also shows a comparison between the ZEUS light curve and 
that derived from the corresponding SPH simulation of the same 
system. The SPH simulation - which we discuss in more detail in 
34.51 - is fully three-dimensional (and hence not optimal for mod- 
elling an axisymmetric problem) but is otherwise set up and anal- 
ysed in a manner that is as close as possible to the ZEUS run. There 
are two main differences. First, individual SPH particles (whose 
smoothing length can vary with time) have a fixed internal energy, 
whereas in the ZEUS runs the sound speed is fixed in Eulerian co- 
ordinates. Second, the SPH viscosity used for computing the im- 
plied dissipation is similar but not identical to that used in ZEUS. 
In particular, it includes both a linear and a quadratic term in the 
divergence of the velocity field (see eq. |24]and eg. I22t. Despite 
these differences we find good agreement between the light curves 
derived from the SPH (dashed lines) and ZEUS (solid lines) sim- 
ulations. The agreement is both in the overall normalization and 
in the rate of decay past the peak. The agreement is not, however, 
perfect. In order to identify the source of the discrepancy we calcu- 
lated the implied SPH dissipation due solely to the quadratic term 
in the artificial viscosity. The resulting light curve is also shown 
as the dotted line in Fig. [8] it is in significantly better agreement 
with the results of the ZEUS runs. This leads us to conclude that 
the different treatment of the artificial viscosity is the main source 
of differences between the two codes predictions. 

We have also tested how sensitive the results are to the magni- 
tude of the artificial viscosity. The potential for unphysical numer- 
ical dissipation requires particularly careful consideration for the 
case of perpendicular kicks, since in this case t he initial perturba- 
tion is an m = warp tha t gives rise to a wave dUubow & Pringlel 
|i993l;lF er reira & Qgilvid FoOS). How such a wave physically dis- 
sipates is not obvious. To test the numerical reliability of the de- 
rived light curves, we have run additional 9 = 90° simulations that 
differ by a factor of two in linear resolution (ZEUS, not plotted) 
or in mass resolution (SPH, shown later in Fig. [T3}. The higher 
resolution runs have lower effective numerical viscosity than the 
lower resolution realizations. For both sets of runs, we find that the 
shape and magnitude of the light curve near the peak is indepen- 
dent of the strength of the numerical viscosity. On the other hand, 
the very early-time behaviour of the ZEUS light curve is found to 
be resolution (and thus viscosity) dependent in the sense expected if 
the dissipation is a numerical phenomenon^ Our conclusion from 
these tests is that the very early rise of the luminosity cannot be 
accurately modeled without the use of a physical model for wave 
dissipation within the disc. The behaviour of the light curve near 
the peak, conversely, appears to be robust and well-captured by 
the code viscosity. We interpret this robustness as being due to the 
much more violent fluid motions associated with the large distor- 
tions seen in Fig.|7] 

4.2 Results for in-plane kicks: razor-thin discs 

In-plane kicks were also simulated using both ZEUS and SPH. It 
is important to realize at the outset that, unlike in the case of per- 

' Higher resolution reduces the luminosity by an amount that is roughly 
consistent with the reduction in artificial viscosity. 



© 0000 RAS, MNRAS 000, 000-000 



8 Rossi et al. 




Figure 7. Visualization of the density in the (r, z) plane (0 < r < 5, -2.5 < z < 2.5) following a kick perpendicular {6 = 90°) to a disc with a surface density 
profile S oc r"'/^ . By the end of the simulation the unbound mateiial has been lost while the remaining bound disc has largely settled back into the equatorial 
(z = 0) plane. The dashed vertical line in each image marks the cylindrical radius r = 1 . 



pendicular kicks, the ZEUS and SPH simulations of in-plane kicks 
model different physical systems. Our ZEUS runs are strictly two- 
dimensional in (r, (p), and thus model a razor-thin disc with a verti- 
cal thickness that is both negligible and constant with radius. Our 
SPH runs, on the other hand, model the effect of the same kick on a 
three-dimensional disc whose scale height is non-zero and increas- 
ing with radius. Although the resulting motions are still predomi- 
nantly confined to the (r, 0) plane, the three-dimensional structure 
of the disc allows gas from large radii to flow inward over the sur- 
face of the disc ballistically without forming a prompt shock as is 
inevitable in a strictly two-dimensional system. As we will show, 
this results in significant changes to the resulting light curve. Here 
we discuss the ZEUS results, deferring the SPH results to 94.51 
where they are presented as part of the investigation of arbitrary 
kick angles. We emphasize that it is not obvious which model is 
closer to physical reality. A real disc is of course three-dimensional, 
as modeled with SPH, but it is also much thinner, at the radii of 
interest, than the disc simulated numerically and hence arguably 
closer to the ZEUS razor-thin limit. 

For the ZEUS simulations of in-plane kicks the computational 
grid extends from r = 0.02 to r = 12, with the initial disc occupying 
the radial range 0.1 < r < 10. For this case, rj, = 0.172 and Tuj, = 
5.8. We use 400 grid points in r, logarithmically spaced to give a 
radial resolution ^ 1.016 at all radii, and 300 grid points in 

(p. All luminosities are scaled to represent a disc mass ratio ^codc = 
6 X 10-*. 

The evolution of the surface density in the kicked disc is 
shown in Fig. |9] for the case of a steeply declining surface den- 
sity profile with p = 3/2. We identify three main phases to the 
evolution, which can be matched to different parts of the bolo- 
metric light curve plotted in Fig. [TO] In the initial phase, a wave 
propagates outward through the disc, depositing energy in annuli at 
successively larger radii. The rate of energy deposition in the sym- 
metric bound region of the disc (r < r^, = 0.172), which includes 
both dissipation of the prompt energy input and accretion energy, 
peaks at t/ty = 0.065. For the p = 3/2 model in which there is 
a substantial amount of mass at small radii, the overall luminos- 
ity peaks shortly afterwards (f/?v = 0.08). The wave then moves 
on through the outer regions of the disc, during which time the lu- 
minosity decays roughly as a power-law with an index i- ^ -0.6. 
This phase ends at about the time when the wave reaches the outer 



edge of the disc. Afterwards, there is an intermediate period dur- 
ing which the total luminosity declines steeply. Finally, there is a 
phase in which the luminosity arises from the infall of low angular 
momentum gas into the inner regions of what remains of the disc. 
Inspection of animations of the simulation shows that the infalling 
gas shocks and enters the disc as dense streams. The inner disc at 
late times is highly non-axisymmetric and time variable. The as- 
sociated dissipation rate is approximately constant but with large 
amplitude fluctuations (see also Fig.ll It. 

The same general behaviour is observed regardless of the 
choice of initial surface density profile, though there are changes to 
the slope of the light curve during the first phase and to the amount 
of low angular momentum material that eventually accretes. Fig.ll II 
shows the energy dissipation rate in a disc that has the same total 
mass but a shallower surface density profile (p = 0.6). As expected, 
the amount of energy dissipated in the outer annuli, whose contri- 
bution peaks at later times, is increased relative to the p = 3/2 disc 
model. As a result, the total luminosity from the disc rises through- 
out the first phase, with a power-law index i- ^ 0. 15, until the wave 
leaves the disc. Since the duration of the energy release is some- 
what more extended the peak luminosity is also reduced. 



4.3 Surface density evolution 

The different behaviour seen for in-plane kicks and perpendicular 
kicks appears to reflect the increased importance of kick-induced 
accretion for small kick angles. This is shown in Fig. [T2l where 
we plot the evolution of the surface density of the bound portion of 
the disc for the in-plane (upper panel) and perpendicular (bottom 
panel) runs. For the perpendicular kick case we find that after the 
kick the inner disc settles back to a surface density profile that is 
very similar to that originally imposed. Note that this settling pro- 
cess takes some time, and it is complete only out to r 0.5 by 
?//v = 10. A negligible amount of gas ends up interior to the origi- 
nal inner edge of the disc. In contrast, there is substantial accretion 
in the = 0° run, and this is evident even at early (t/ty =0.1) times 
when the light curve is undergoing its initial rise. 
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Figure 9. Rendering showing the evolution of the surface density of the disc following an in-plane kick. In this high resolution simulation the disc extends 
up to r = Tub. During the peak phase of energy dissipation (left-hand panel) energy is dissipated at successively larger disc radii as an outward moving wave 
propagates through the gas. The outer part of the disc is unbound and escapes ballistically. After the wave reaches the outer edge of the disc (centre panel) 
the rate of decay of the energy dissipation rate steepens markedly. At late times (right-hand panel, spatial and colour scale adjusted to show structure in the 
innermost regions) low angular' momentum gas continues to accrete on to the bound remnant of the original disc, releasing energy at a low level and forming 
a highly non-axisymmetric accretion flow. 



10-1 



10-5 



-o 



e 



10- 



10- 




10-1 



10-5 



10- 



10- 




Figure 10. The implied energy dissipation rate as a function of time for an 
in-plane kick into a razor-thin disc with a surface density profile Z oc r"''^ 
and a fixed sound speed profile c, oc r"^!'' . The disc extends up to r = 10. 
Time is measured in units of the dynamical time (fy ) at , while the energy 
dissipation rate is in code units for a disc mass Md/M = 6 X 10"* between 
r = 0.1 and r = 10. The solid black curve shows the total energy dissipation 
rate evaluated across the whole computational domain, while the dashed 
black line is the dissipation rate for 0. 1 < r < 10 for comparison with the 
SPH results. This latter has been divided up into individual contributions 
from annuli at successively greater radii: 0.1 < r < 0.172 = rj, (red short- 
dashed line), 0.172 < r < 0.55 (blue long-dashed line), 0.55 < r < 1.79 
(green dot-dashed line) and 1.79 < r < 5.8 = r„5 (cyan dot-long-dashed 
line). The dotted black curve (at the bottom) shows the result for a disc with 
no kick but an instantaneous black hole mass loss of 0.03M. 



Figure 11. The implied energy dissipation rate as a function of time for a 
disc with a flatter surface density profile S oc r"^'^ . The total energy dissi- 
pation rate (shown as the solid black curve) peaks at a smaller value than 
for S oc r-'^'^. Light curves for individual annuli are plotted as in Fig. 1101 
except that no mass-loss light curve is shown here. 



4.4 Mass-loss simulation 

Fig. [To] also shows how the energy dissipation rate in the disc due 
to the kick compares to that which would arise as a result of the 
black hole mass loss experienced upon merger (the dotted line at 
the bottom of the figure). For simplicity we set up a disc model 
identical to that used for the kick simulation, but set V = and 
imposed an instantaneous mass loss of 3%. The simulation was run 
until ijtM = 2. Given our parameters, the energy dissipation rate in 
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Figure 13. Light curve obtained for 6 = 30° and 6 = 90°, from SPH 
simulations at two different resolutions: 4 million particle run (red dashed 
lines) and 8 million particle run (solid black lines). 



Figure 12. Evolution of the disc surface density with time from ZEUS sim- 
ulations of kicks at 6 = 0° (upper panel) and = 90° (lower panel) to 
the disc plane. For clarity only gas that remains bound to the black hole is 
included in the calculation of the surface density. Both discs had initial sur- 
face density profiles Z oc r"^'^. (solid lines, plotted here in arbitrary units). 
Additional curves show the profile at t/t\ =0.1 (dotted lines), t/t\ = 1 
(short-dashed lines) and at the end of the simulations (long-dashed lines 
plotted at f/fy = 20 for the in-plane run and t/t\ = 10 for the perpendicular 
kick). 

the disc induced by the mass loss peaks early (t/t\ = 0.085) at a 
luminosity that is more than two orders of magnitude below that 
generated by the kick. We can therefore consistently ignore black 
hole mass loss in our simulations. Of course this does not mean that 
mass loss is always physically negligible, since (to quote a trivial 
example) mass loss must occur in mergers whose symmetry is such 
that no recoil is produced. The analytic discussion in Section 1231 
provides an estimate of when mass loss might be competitive as a 
source of disc energy dissipation. 

4.5 Results for out-of-plane kicks 

In this Section, we discuss the results of SPH simulations of discs 
subject to kicks at arbitrary angles. Our simulations are for kick 
angles 9 = 0° (with a finite thickness disc), 6 = 15°, 8 = 30°, 6 = 
60° and 6 = 90° . As previously, we adopt a nominal disc to black 
hole mass ratio q = 6x 10 As an initial numerical check we show 
in Fig.[T3]the light curves for 9 = 30° and 9 = 90° at two different 
resolutions: 4 million particles (dashed line) and 8 million particles 
(solid line). Excellent agreement is obtained, providing evidence 
that our simulations have converged. This result, together with the 
results of the code comparison ai9 = 90° , makes us confident that 
any uncertainties in the numerical computation of the light curves 
from the SPH simulations are small. For practical purposes, they 
are certainly much smaller than physical uncertainties arising, for 
example, from poor knowledge of the actual initial surface density 
profile at the relevant radii. 

Visualisations of the evolution of the column density from the 
9 = 30° run are shown in Fig. [14] (edge-on view) and in Fig. [15] 



(face-on view). For this kick angle the two characteristic radii are 
approximately equal to x 0.2 and r^b ~ 4.8, respectively. Al- 
though the geometry is now more complex the overall behaviour 
resembles that seen in the special cases discussed earlier. As the 
black hole is ejected from the centre of the disc it carries with it 
the innermost part of the disc, while the outer part lags behind and 
is dispersed. The bound mass in the simulation agrees with that 
predicted analytically and is constant with time. The edge-on im- 
ages show that there is a brief initial phase during which the disc is 
warped, but as was found for 9 = 90° ,the gas quickly readjusts to 
orbit within a single plane. In this case, however, the plane is tilted 
with respect to that of the original disc (see discussion in 32.2| l. 
The face-on images show that the final extent of the remnant disc is 
comparable with r\,. As we will see in the following section, most 
of the dissipation occurs in this region (r rj,). 

The implied rate of energy dissipation as a function of time is 
plotted in Fig. [T6] with a black solid line. The contribution to the 
light curve from gas at dilferent distances from the black hole is 
shown in the same figure. We have divided the disc into three broad 
annuli: 0.1 < r < r^, ry, < r < rub and r > rub. Within each an- 
nulus, we computed the energy dissipation rate by summing over 
all particles that are instantaneously within that region. As was the 
case for the in-plane kick, we see an initial rise in the luminosity 
as the compression wave propagates outwards and energy is dis- 
sipated at successively larger radii. The peak in the total curve at 
t ^ 0.9 is primarily due to dissipation in the inner region (r < rb), 
although material in the intermediate zone rb < r < rub - whose 
contribution dominates the light curve during its declining phase - 
is responsible for a substantial amount of energy when integrated 
over time. Within this ring (rb < r < r^b) most of the dissipation 
occurs close to rb. The contribution to the light curve from gas at 
r > rub is negligible and consistent with zero, as we assumed in the 
analytic calculatior[3. 

We also note that both the total light curve and the light curves 
of the two individual rings within rub show some evidence for a 
two-component structure. The presence of two components is more 

* We verified in particular that particles originally at r > Tub are lost with- 
out encountering shocks or compressions that result in significant heating. 
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Figure 14. Visualization of the column density (in code units) at different times, from a simulation with a kick at = 30°. Part of the disc remains bound to 
the black hole, and this material rapidly settles into a relatively stable tilted configuration. Unbound gas can be seen leaving the system in the upper portion of 
the panels. 




clearly seen in the d = 15° light curve, shown as Fig. [T7] For this 
figure we have computed the dissipation not only as a function of 
the current Tudius of particles within the disc (upper panel), but also 
by binning the particles depending upon their initial location (lower 
panel). This allows us to (at least partially) separate the contribu- 
tions from the prompt deposition of kick energy and that arising due 
to accretion. We find that most of the dissipation occurs within r^ 
(upper panel). At early times this dissipation is partially associated 
with gas that originated in this region (lower panel), and we thus 
infer that the initial rise (up to and through the first shallow peak) 
in luminosity is due to both prompt dissipation of kick energy and 
accretion that involves small changes in orbital radius. The second 
more prominent peak, on the other hand, occurs as energy is dissi- 
pated within Tb (upper panel) in gas that was mostly originally in 
the region between rj, and r^i, (lower panel). We identify this con- 
tribution as being almost entirely due to accretion of low angular 
momentum matter that falls from larger radii and deposits energy 



at smaller radii. The deposition occurs via mild shocks and this may 
cause matter to lose additional angular momentum, sink further to- 
ward the centre and therefore release additional potential energy. 



In order to quantify whether accretion occurs beyond the level 
expected from angular momentum conservation, we have made use 
of the Lagrangian nature of SPH: we plot the final (t/ty = 8.5) ra- 
dial position of bound particles as a function their initial position on 
a particle-by-particle basis (Fig.|18ll. The data are from the 9 = 30° 
simulation. In Fig.[T8] the red line shows the locus where particles 
should lie if they conserved their angular momentum (equation |7]l. 
We find that to a first approximation angular momentum conserva- 
tion is a good assumption. However, there is clearly a trend for the 
bound gas - especially from rj, < r < r^b - to fall deeper into the 
potential well, and form a remnant disc that is smaller than would 
be predicted under our assumption. This is the result of angular 
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Figure 16. Light curve derived from an SPH simulation of the 6 = 30° 
case. dE/dt is in code units. Here we also show the individual contributions 
arising from energy dissipated in thi'ee different radial regions as marked. 
Clearly most of the dissipation occurs within 
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Figure 18. Final (at t/t\f = 8.5) radial position of a bound (e < 0) SPH 
particle in units of rb as a function of its initial radial position in the same 
units, for the d = 30° case (blue symbols). For clarity, only 10"' of the total 
bound particles are plotted. For comparison, we plot as the red soUd line our 
analytic expectation (equation|2). The black solid line marks where r = r^- 
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Figure 17. Light curve for 8 = 15° in code units. This particular simulation 
used 4 million particles. Upper panel: The total dissipation rate (black soHd 
line) is the sum of the energy dissipated within ri, (red short-dashed line) 
and between rj, and r^), (blue long-dashed line). These contributions are 
calculated summing over particles present at that time in the region. Lower 
panel: the same as above, but here we show the individual contributions 
of those particles initially (at 1 = 0) within rb (red short-dashed line) and 
between f-b and f-yb (blue long-dashed line). 



momentum mixing, where particles with larger initial angular mo- 
mentum transfer part of it to particles with an initially lower one. 

The final formation of a relatively compact disc can also be 
seen in Fig. [19] where we plot the evolution of the surface density 
for this run. Most of the bound particles in the final disc end up 
within Tb, with a steep tail at larger radii that drops to almost zero 
by about 3rb. This behaviour does not match that predicted from 
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Figure 19. Surface density (code units) as a function of radius. At the be- 
ginning of the simulation all particles are distributed with a r"3/2 pj-Qfjig^ 
while the bound particles follows a steeper profile of Z oc r"". At the end 
of the simulation {tlty = 8.5) the bound particles are mostly within Sb ™d 
their surface density distribution falls-off rapidly for larger radii. 



the model assuming circularization at constant angular momentum. 
Rather, there is an indication of angular momentum transfer and ad- 
ditional accretion of matter (see also Fig.llSt. Since exact conser- 
vation of particle angular momentum was assumed in our analytic 
model, these results imply that it cannot provide a precise quanti- 
tative estimate of the total energy release, but rather (for this kick 
angle) a lower limit. Longer duration simulations will be needed 
to determine how much energy in excess of our simple model is 
indeed dissipated. 

The dependence of the shape and amplitude of the light curve 
on the direction of the black hole recoil is illustrated in Fig. [20] The 
peak luminosity and the total dissipated energy increase for kicks 
that are more closely directed toward the disc plane. This trend is 
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evident for angles of 30° and smaller (the Q = 90° and 6 = 60° runs 
are quite similar). The shape of the light curve is a complex function 
of the kick angle. We interpret the behaviour seen in the simulations 
as being a consequence of the existence of two components to the 
light curve, whose separation in time is a function of kick angle. 
One contribution is associated with gas in the inner disc (within 
and close to r^,). As the kick angle decreases r^, also shrinks, and 
the energy dissipated within this region of the disc is released on a 
shorter time scale. A second contribution, due to accretion, occurs 
on a time scale that is related to r^i,, which increases for kicks that 
are directed closer to the disc plane. The superposition of the two 
components results in a light curve that is clearly double-peaked 
only for a narrow range of kick angles, but whose evolution with 
kick angle qualitatively matches that seen in the simulations. 



We also note that the result for the 6 = 0° simulation presented 
here, for a finite thickness disc, difi'ers significantly from that ob- 
tained with ZEUS for a razor-thin disc 34.21 (compare Fig.l20lwith 
Fig.llOb. In particular, for the finite thickness disc we do not observe 
the very rapid rise to peak luminosity seen in the ZEUS simulation. 
Instead, the 9 = 0° light curve continues the trend toward longer 
time scale energy release that is seen in the 9 = 30° and 9 = 15° 
simulations. Analysis of the SPH simulation suggests that the rea- 
son for the difference between the two results lies in the presence 
of three-dimensional flows of gas across the surface of the kicked 
disc. They cannot occur in a strictly two-dimensional calculation. 
In the SPH runs the kick causes some of the gas in the outer disc 
(at relatively large height above the mid-plane) to flow inward over 
the surface of the disc, before releasing energy at smaller radii and 
later times. Inspection of the vertical profile of energy dissipation 
in the SPH run shows a substantial release of energy in the surface 
layers of the disc. 



Finally, we compare in Fig. [21] the cumulative energy release 
as a function of with the analytic estimate discussed in Sec- 
tion |2.5l The figure shows (as the red stars) the total energy release 
derived by integrating the high and moderate 6 light curves shown 
in Fig.|20]through to the end of the simulations. It is clear that en- 
ergy dissipation is still ongoing at the (arbitrary) epoch at which 
the simulations stop, so our estimates of the total release are lower 
limits. Nonetheless, the simulations clearly show a trend toward 
increasing energy release for smaller kick angles that - for kick 
angles between 90° and 15° - agrees surprisingly well with that 
predicted analytically. For the 0° run (whose integrated energy re- 
lease is not plotted because the simulation was of shorter duration) 
the light curve (Fig. 1201) shows a further increase in total energy as 
compared to the 15° run, though we do not find evidence for the 
very large (in fact divergent) increase predicted analytically. We at- 
tribute this difference as being due to the fact that the analytic diver- 
gence arises because a small amount of mass in the disc is predicted 
to circularize at zero radius. Any hydrodynamic mixing suffices to 
give the initially zero angular momentum gas a small amount of 
angular momentum, eliminating the divergent behaviour. This re- 
sults in a smaller total energy release. We caution, however, that 
the maximum value of the total energy release at small kick angles 
may well depend upon the physical thickness of the disc, which is 
larger in the simulations than in real systems. For very thin discs 
the analytic trend toward increasing energy release might extend to 
smaller kick angles than those obtained here. 
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Figure 20. Total (from the whole disc) rate of energy dissipated as a func- 
tion of time for different recoil directions. The figure shows results from 
SPH simulations at 6 = 90°, 9 = 60°, 8 = 30°, = 15°, together with the 
in plane case (6 = 0°). All curves have the same disc model with a surface 
density profile Z oc r^^l^ and a total mass (relative to the black hole) of 
q = 6x 10"'' between r = 0. 1 and r = 10. 
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Figure 21. Comparison of the analytic prediction for the total energy release 
(solid curve) with the simulation results (stars), as a function of kick angle. 

5 DISCUSSION 

Our numerical results can be scaled to represent a wide range of 
different systems. In doing so, three parameters enter the problem: 
(1) the black hole mass, which affects the characteristic time scale 
tv over which energy is deposited into the disc, (2) the kick ve- 
locity, which affects both the time scale and the normalization of 
the luminosity, and (3) the disc mass, which alters the luminos- 
ity. The effect of kicks on surrounding discs might in principle be 
observable in two distinct regimes. For low-mass black holes with 
M ~ 10^ Mq the time scales for the initial rise in the light curve in 
the source frame range from ~ 10 yr for V = 300 km s"' down to as 
little as ~ 1 yr for V > 10^ km s"'. It may therefore be possible to 
detect a time-variable electromagnetic counterpart from the kicked 
disc following a merger that has been approximately localized by 
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gravitational-wave (or other) observations ( iLippai et alj|2008h . For 
high-mass black holes (M ~ lO^Afg) the long time scales for vari- 
ability and the rarity of mergers in this mass range combine to make 
such triggered searches unfeasible. It has been suggested, however, 
that the spectrum of a kicked disc might be sufficiently distinct from 
that of other sources as to allow identification in wide-area surveys 
(Schnittman & Krolik 2008). 

In addition to the black hole mass and kick velocity, the disc 
mass is an important parameter, since the luminosity is propor- 
tional to the surface density at the characteristic radius Ry. It is 
important to note that this is a large radius: for M = 10* Mq and 

V = 300 kms"', we have Ry = 0.05 pc (10* GM/c~), while 
for M = 10'* Mo and V = 10' km s"' we have Ry = 0.43 pc 
(» 10^ GM/c^). At these radii an argument can be made that the 
maximum surface density of a geometrically thin accretion disc is 
limited by the onset of fragmentation due to the disc's self-gravity 
(Gamm ie 2001; Goodman 2003; Rice, Lodato & Armitage 200^ 
lRafikovH2005l ; Ifor a review see Lodatol |2007|) . A simple estimate 
can be derived by m aking use of the steady-state models calcu- 
lated by iLevinI ilOOt . In these models the maximum value of Z 
is a function only of fl, and varies from E^ax ~ 2 x 10' g cm"^ 
at a radius where the orbital period P = 10^ yr down to E„,ax ~ 
40 g cm"^ or less at f = 10' yr. Writing the maximum disc mass 
as M„ax ~ ^Ry^il^v), where l.(Rv) is the maximum value of the 
surface density that would be stable against fragmentation at Ry, 
we find that for a range of kick velocities between 300 km s"' and 
10' km s"' the disc around a 10* Mq black hole would be unlikely 
to exceed q ~ 10"' of the black hole mass. In scaling our results to 
M = 10* Mo, we therefore use the value 17 = 6x 10 discussed ear- 
lier. For a 10^ Mo black hole a somewhat more massive disc could 
be stable, and for this case we scale to q = 6x 10"'. 

Given the important role that the assumed disc mass plays in 
determining the observability of emission from kicked discs, it is 
worth stressing that these estimates are crude. There are few ob- 
servational constraints on discs at sub-pc scales in galactic nuclei. 
Theoretical work only strictly excludes the existence of very mas- 
sive discs that are heated by viscous processes and cool radiatively. 
Therefore, it would not surprise us if our "limits" could be ex- 
ceeded. 

Fig. |22] shows examples of the scaled light curves for three 
sets of parameters: (1) a baseline case with M = 10* Mq, V = 
300 km s"' and q = 6 x lO"''; (2) a more optimistic (as far as 
observability goes) case with M = 10* Mo, V = 10' kms"' and 
q = 6 X 10"**; and (3) a high mass case with M = lO'^Mo, 

V = 10' km s"' and q = 6x 10"'. These parameters are similar to 
those considered by Schnittman & Krolik (2008) except for the fact 
that, for the reasons outlined above, we use a substantially smaller 
disc mass. For each parameter set we plot the total luminosity from 
the 6 = 15° and 9 = 90° simulations, which roughly bracket the 
range of likely behaviour seen in our simulations. 

The maximum luminosity obtained from the M = 10* Mq 
models is L 10*' erg s"', which corresponds to about 10% of 
the Eddington luminosity Ledd = 1-3 x lO** erg s"' for this mass of 
black hole. Reaching this luminosity, however, requires a combina- 
tion of circumstances that may be uncommon: a high velocity kick 
{V = 10' km s"'), directed almost into the plane of a disc, with a 
steep surface density profile. The peak luminosity for most of the 
other models plotted in Fig. [20] is smaller: around 2 x 10*^ erg s"' 
(1.5% Z-Edd) for V = 10' km s"' and 6 - 90°, and approximately 
lO'"' ergs"' (lO-^LEdd ) for the models with V = 300 kms"'. 
These values encompass the ra nge of luminosities (6.3 x lO^'^Lgdd 
to 1.6 X lO^^Lsdd) estimated bv lLippai et alj ( I2OO8I) . 
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Figure 22. The predicted bolometric light curve of the kicked disc after 
scaling the numerical results to represent three different classes of systems. 
The solid curves show the light curves for 6 = 15° (upper curve) and d = 
90° (lower curve) runs, following scaling to M = 10*A/o, V = lO' km s"', 
and a disc mass ratio g = 6 X 10"''. Both runs are for a surface density 
profile with p = 3/2. The short dashed curves are for identical parameters 
except for a lower kick velocity of V = 300 km s"' . The long dashed curves 
are for a system with M = 10* Mq, V = lO' km s"', and a disc mass ratio 
q = 6x 10"'. 



Fig. |23] shows the predicted bolometric flux from the disc as 
a function of source redshift z- We assume a luminosity distances 
appropriate for a standard cosmology. As expected from the fact 
that even the most luminous sources are sub-Eddington, none of 
our sources are predicted to be very bright. Electromagnetic coun- 
terparts might be detectable if V lies toward the upper end of the 
considered range, especially if the merger is nearby. Low veloc- 
ity kicks (V = 300 km s"'), on the other hand, yield very low lu- 
minosities that would not be detectable at plausible cosmological 
distances. 

Fig.|22]also shows the predicted light curve for a disc around a 
10** Mo black hole following a kick of 10' km s"' . These parameters 
have been chosen to match those adopted by Schnittman & Krolik 
(2008), although our curves are calculated for a steeper surface den- 
sity profile than the S oc r"''' that they used. Nonetheless we find 
that (for the 9 = 15° run) the time of the predicted peak (? ~ 10' yr 
in the source frame) is similar to that derived by Schnittman & Kro- 
lik (2008). The amplitude, however, is more than an order of mag- 
nitude lower (about lO*"* erg s"' as opposed to a few 10"*^ erg s"'). 
This difference arises because we have limited the disc mass to the 
maximum that is allowable before self-gravity would result in frag- 
mentation. Unless there is some way to circumvent this limit, we 
conclude that the prospects for detecting the emission from kicked 
discs in surveys are poor. 

In addition to the energetic considerations discussed here, the 
actual observability of emission froin kicked discs will also depend 
upon the spectral band in which the radiation is emitted. We expect 
emission in the soft X-ray, if the disc is optically thin at the radii 
of interest, or if the bulk of the energy is ultimately deposited in 
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Figure 23. The predicted bolometric flux from the disc as a function of the 
redshift z of the source is plotted for two kick velocities, V = 300 km s" ' 
and V = km s"'. For each velocity we show two curves that approx- 
imately encompass the range of peak luminosities implied by our simula- 
tions. The lower dashed curves show results for typical kick orientations 
(corresponding to luminosities in code units of 10"'), while the upper solid 
curves show results for more optimistic geometries in which the kick is 
directed close to the disc plane (IQ— * in code units). 

the upper layers of the disc atmosphere The emission would be at 
a temperature of Tx ~ 1 keV, that correspond s to the post-shock 
value behind a shock of a few hundred km s ' jLi ppai et al.ll2008l : 
[Shields & Bonning|l2008l) . Energy deposition at the midplane of an 
optically thick disc, conversely, would result in thermal emission 
in the infrared (Schnittman & Krolik 2008). The results of iLevinI 
( l2007l) imply that self-gravitating discs ought to be optically thick 
out to radii where the orbital period is P ^ 600 yr. Hence the high- 
velocity kicks that are most observable would deposit their energy 
at radii where the disc was optically thick. We caution, however, 
that this does not necessarily imply that all of the energy is ther- 
malized into the infrared. Waves within discs can - under certain 
circumstances - deposit energy preferentia lly in the low den sity at- 
mosphere rather than at the disc midplane l lBate et aTll 20021) . More 
detailed simulations, that include a realistic treatment of the disc's 
vertical structure, will be needed to determine where energy is de- 
posited, and what the resulting spectral signatures are. 



6 CONCLUSIONS 

In this paper we have investigated the effect of a post-merger kick 
on the dynamics of a geometrically thin accretion disc, surrounding 
the newly merged black hole. The existence and strength of kicks 
following black hole mergers are now secure predictions of General 
Relativity. Whether small-scale gas discs typically attend mergers 
remains, however, uncertain. If gas discs are commonly present the 
perturbation to the gas caused by the kick will generate an electro- 
magnetic counterpart to the merger, which may be detectable if its 
amplitude and time scale are observable at cosmological distances. 



The main conclusions of our study are: 

1. The dimensional assumption that the magnitude of energy re- 
lease is (l/2)K'Sv^v ™ underestimate. The energy available for 
dissipation varies strongly with the angle between the kick and the 
disc plane, and for kicks close to the plane our analytic estimate 
exceeds (l/2)K^Ev./?v W up to three orders of magnitude. A large 
increase in the energy release for kicks close to the plane is con- 
firmed numerically. 

2. For most orientations of the kick, accretion energy (i.e. energy 
liberated when gas in the disc loses angular momentum and falls in- 
ward) dominates over the direct energy input to the gas in the frame 
of the kicked black hole. Accretion energy is particularly important 
for kicks that are directed toward the disc plane. The importance 
of accretion can be demonstrated analytically, and is confirmed by 
numerical simulations. 

3. We have run SPH and (for the special cases that are two- 
dimensional) finite difference numerical simulations to investigate 
the evolution of the disc and the rate of energy dissipation, follow- 
ing a kick at an arbitrary angle to the disc plane. The simulations 
yield explicit predictions for the form of the light curves as a func- 
tion of both kick angle and surface density structure within the disc. 

4. The observability of emission from kicked discs depends upon 
the kick velocity, the orientation of the kick relative to the disc 
plane, and the mass of the disc at the radii where the energy is 
deposited. We have argued that the mass of the disc is limited by 
the requirement that the disc remains stable to fragmentation due to 
self-gravity. If this is true, the decreased disc mass (as compared to 
that assumed in prior works) outweighs the effect of the increased 
energy release per unit mass. As a result, the predicted luminosity 
around massive black holes is substantially smaller than previously 
thought. It is unlikely that such sources can be identified via wide- 
area sky surveys. 

5. The most feasible observational probe of the phenomena dis- 
cussed here is via identification of variable disc emission follow- 
ing black hole mergers detected by other means (e.g. via detection 
of gravitational waves). Disc emission may be detectable provided 
that the kick velocity is large, especially in the case where the kick 
is directed close to the disc plane. To quantify, a kick velocity of 
10' km s"' offers promising possibilities, but 300 km s"' appears 
to be very difficult or impossible to detect. 

Physically, the magnitude and orientation of the kick is deter- 
mined by the magnitude and direction of the spin of the black 
holes immediately prior to their merger. Our results emphasize 
that the observability of this class of electromagnetic counter- 
parts depends sensitively on the distribution of the pre-merger 
spins. Studies, such as that bv Bogdan ovic et al] ( |2007|) and that by 
|King, Pringle & Hofmann (2008), that attempt to predict the evolu- 
tion of black hole spin during the earlier phases of merger are thus 
particularly important for determining whether counterparts will be 
detectable. Determining the spectral signature of the kicked disc is 
also important, and this will require simulations of the disc evolu- 
tion that include more complete treatments of the disc physics. 
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